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The Role of Gas In Maintaining Quasi-Steady Spiral Structure in 
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ABSTRACT 

We study the dynamical evolution of spiral structure in the stellar disks of isolated galaxies 
using high resolution Smoothed Particle Hydrodynamics (SPH) simulations that treat the evo- 
lution of gas, stars, and dark matter self-consistently. We focus this study on the question of 
self-excited spiral structure in the stellar disk and investigate the dynamical coupling between 
the cold, dissipative gaseous component and the stellar component. We find that angular mo- 
mentum transport from the gas to the stars inside of corotation leads to a roughly time-steady 
spiral structure in the stellar disk. To make this point clear, we contrast these results with oth- 
erwise identical simulations that do not include a cold gaseous component that is able to cool 
radiatively and dissipate energy, and find that spiral structure, when it is incipient, dies out 
more rapidly in simulations that do not include gas. We also employ a standard star formation 
prescription to convert gas into stars and find that our results hold for typical gas consumption 
time scales that are in accord with the Kennicutt-Schmidt relation. We therefore attribute the 
long-lived roughly time steady spiral structure in the stellar disk to the dynamical coupling 
between the gas and the stars and the resultant torques that the self-gravitating gaseous disk 
is able to exert on the stars due to an azimuthal phase shift between the collisionless and 
dissipative components. 
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1 INTRODUCTION 

Spiral galaxies present two faces to observers. When seen in blue 
light, these galaxies look disordered, with spurs, branches and 
feather-like structures seen even in grand design galaxies, while 
K-band images show that these galaxies possess smooth, usually 
two-armed spiral patterns. We study here the dynamical interac- 
tion between the Population I component seen in blue light and the 
older Population II component revealed in the longer wavelengths, 
specifically with the intent of understanding how the interaction 
between the gas and the stars affects the time evolution of spiral 
structure in the stellar disk. 

The idea that the stellar spiral structure is not a material phe- 
nomenon, but a roughly time-steady density wave rotating at a pat- 
tern speed was developed in the seminal work of Lin & Shu (1964; 
henceforth LS64). Toomre (1964) gave the minimum radial stellar 
velocity disperion required to suppress the onset of local instabili- 
ties in a razor-thin stellar disk stabilized by rotation and the velocity 
dispersion of stars on large and small scales respectively. Much of 
the early numerical work of the study of spiral structure utilized 
N-body simulations to study the growth of spiral instabilities in the 
collionsionless components of galaxies. Transient spiral structure 
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was evinced in the early N-body simulations of Hohl (1970), as 
well as in the later works of Sellwood & Caiiberg (1984; henceforth 
SC84). Toomre (1981) showed that a shearing disk with a supply 
of leading waves will swing amplify these leading disturbances to 
trailing disturbances, and was one of the earliest works to empha- 
size the role of cosmological accretion of gas. Tidal disturbances or 
minor mergers can torque the stellar disk and create transient spiral 
structure (Hernquist & Mihos 1995). Some of these early influential 
works noted that the dissipative nature of gas may provide a self- 
regulatory influence in the combined gas-disk system (Lin & Shu 
1964; Bertin & Lin 1996), with dissipation in the gas layer serving 
to continually regenerate the shearing spiral instability (Goldreich 
& Lynden-Bell 1965). 

Contemporary numerical investigations of spiral structure 
have focused mainly on the response of a gas disk when driven by 
a spiral potential, to varying levels of physical detail and numerical 
sophistication. Kim & Ostriker (2002) carried out local magnetohy- 
drodynamic (MHD) simulations to study the development of spiral 
substructure in a gas disk forced by a spiral potential. Chakrabarti, 
Laughlin & Shu (2003) (henceforth CLS) studied the development 
of spiral substructure due to ultraharmonic resonances using global 
simulations which incorporate the self-gravity of the gas and found 
that inter-arm features of long azimuthal extent occurred at reso- 
nance locations, which they called "branches". CLS also found that 
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when the gas disk is forced by a time-steady, but nonlinear spiral 
potential (which is composed of many Fourier components), over- 
lap of resonances produced a highly time-dependent response in 
the gas which they speculated may signal the onset of turbulence in 
the purely hydrodynamic regime. Shetty & Ostriker (2006; hence- 
forth SO06) carried out global MHD simulations and found that the 
magneto-Jeans instability is effective at forming spiral-substructure 
in low-shear (i.e., near the spiral arms) regions when the magnetic 
field is large. The stabilizing effect of magnetic pressure allowed 
these simulations to be carried out for longer times than those of 
CLS. Yanez et al.'s (2007) hydrodynamic simulations also evince 
the role of ultraharmonic resonances in the development of spiral 
sub-structure in the gaseous component. A significant work on the 
analysis of spiral structure in self-gravitating disks is that of Laugh- 
lin et al. (1998) - these authors rigorously show that the primary 
unstable mode is endemic to the disk and not due to the numerical 
implementation. 

One of the essential points evinced by these simulations of 
the response of a self-gravitating gas disk to a time-steady stellar 
potential is that low Q gas disks display a very time-dependent re- 
sponse and continue to do so over many dynamical times (where 
we refer to the Toomre Q parameter, i.e., Q — Cin/trGT,, where a 
is the sound speed of the relevant component; Toomre 1964). This 
is in marked contrast to purely stellar disks which may exhibit tran- 
sient spiral structure (if they are low Q disks) but the spiral patterns 
disappear after a few dynamical times. 

Block et al. (1994; 1996) found that many spiral galaxies have 
a grand design structure in infrared light (Elmegreen & Elmegreen 
1984; Rix & Zaritsky 1995) while appearing quite fiocculent in blue 
light. They interpret this discrepancy in morphology to mean that 
the Population I and Population II components are dynamically de- 
coupled. Recent multi-wavelength observations have begun to ad- 
dress in detail the long term evolution of spiral structure in nearby 
galaxies. Notably, Kendall et al. (2008) infer from analysis of high 
spatial resolution multi-wavelength photometry of M8 1 (from the 
B-band to IRAC bands) a phase shift between the gas and the stars 
that is suggestive of a long-lived spiral structure in the stellar disk. 
Global torques may have a role in fueling star formation over many 
Gyrs if they are able to transport the large amounts of mass in ex- 
tended HI disks (Wong & Blitz 2002) to the inner regions. These 
issues require a global analysis that treats the interplay between gas 
and stars rather than treating the response of the gas to an imposed 
stellar potential as has been studied previously (CLS; SO06, among 
other previous papers). 

In this paper, we employ high resolution three-dimensional 
simulations performed with the Smoothed Particle Hydrodynamics 
(SPH) code GADGET- 2 (Springel 2005) to simulate the evolution 
of self — excited spiral structure in the stellar disks of galaxies 
that are evolving in relative isolation. The recent analysis of the 
chemical and star formation history of the Milky Way by Colavitti 
et al. (2008) suggests that the Milky Way has not undergone recent 
(i.e., after z ~ 1) significant merger events. We study here the prob- 
lem of the generation and evolution of spiral structure in an isolated 
galaxy with gas, stars, and a dark matter halo, where these compo- 
nents are free to interact dynamically with each other. Specifically, 
by including not only the stars with their gravitational forces and 
velocity dispersions, but also the cold interstellar gas with its self- 
gravity, pressure, its ability to cool radiatively and dissipate energy 
through shocks, and form new stars, we focus our study on the dy- 
namical coupling between the gas and the stars and the torques that 
the gas is able to exert on the stars, an effect that has not been 
analyzed in previous studies which have focused on the effect of 



the stars (whether live or modeled through a potential) on the gas 
(Hernquist & Mihos 1995; CLS; S006; Kim & Ostriker 2007). We 
also comment on the coupling with the live halo and its role as a 
sink of angular momentum, as has been noted in previous studies 
of bar instabilities (Athanassoula 2003; Martinez- Valpuesta et al. 
2006; Berentzen et al. 2007). We find that the dissipative nature of 
the gas leads to sharp, asymmetric azimuthal features, with the 
gas response leading the stellar response, which allows for a phase 
shift between these components such that the gas is able torque the 
stars. It is this angular momentum transport interior to co-rotation 
that allows for the stellar disk to be unstable to spiral structure for 
a longer time (many revolutions) compared to simulations that do 
not include the dissipative gaseous component. 

The paper is organized as follows. In §2, we review the sim- 
ulation methodology, specifically the initial equilibrium structure 
of the disks and implementation of radiative cooling, viscosity, and 
star formation prescriptions for turning gas into stars. In §3, we give 
our results for simulations that treat the collionsionless components 
only, and in §4, we describe results from simulations that include 
the time evolution of gas, stars, and a dark matter halo. We discuss 
our results in the context of observations in §5, and we conclude in 
§6. The Appendices are devoted to a detailed exploration of simu- 
lation parameters and assumptions that can affect the gas-star an- 
gular momentum transport, specifically the artificial bulk viscosity 
and equation of state for the gas. 



2 SIMULATION METHODOLOGY 

We employ the parallel TreeSPH code GADGET-2 (Springel 2005) 
to perform simulations of isolated disk galaxies. GADGET-2 uses 
an N-body method to follow the evolution of the collionsionless 
components, and SPH to follow the gaseous component. Energy 
and entropy are both conserved in regions free of dissipation. The 
simulations reported here (unless otherwise noted) have gravita- 
tional softening lengths of 25 pc for the gas and stars, and 50 pc 
for the halo. The number of gas, stellar, and halo particles are 10 5 , 
1.08 x 10 6 and 10 6 respectively for our fiducial case. The total an- 
gular momentum of these simulations is conserved to 0. 1 % over 3 
Gyr for a time integration error tolerance accuracy of 0.025. In Ap- 
pendix A, we give some results of a few cases where we vary the 
artificial bulk viscosity and the equation of state - these simulations 
are performed at lower resolution and their specifications are noted 
in the appendix. 

The implementation of viscosity in GADGET-2 is described 
in detail in Springel et al. (2001) and Springel (2005), so we briefly 
review the results here. GADGET-2 uses a bulk artificial viscos- 
ity which is the shear-reduced version (Balsara 1995; Steinmetz 
1996) of the standard Monaghan and Gingold (1983) artificial vis- 
cosity. Specifically, the form of the artificial viscosity experienced 
between particles i and j is taken to be given by: 

a (cj + cj - 3wj) Wj 
2 pij 

where pij is the arithmetic mean of the SPH density estimator 
(equation 5 in Springel 2005) of particles i and j, a and Cj are the 
sound speeds, and w t and Wj are the relative projected velocities. 
Here, we adopt a fiducial value of 0.75 for the bulk artificial viscos- 
ity a parameter. This value has been empirically recommended by 
Springel et al. (2001) by comparing the results of the shock tube test 
to other numerical implementations. In discontinuities that develop 
in the gaseous flow, the entropy generation and width of the shock is 
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proportional to this parameter. We vary this parameter in Appendix 
A and find that angular momentum transport from the gas to the 
stars is inhibited for very low values of the artificial bulk viscos- 
ity parameter. For low a values, the effect of dissipation in shocks 
is prohibited such that the gaseous component becomes more akin 
to the collionsionless component and the stars no longer exhibit a 
significant time lag or azimuthal phase shift in their response. 

Radiative cooling and heating are implemented in GADGET-2 
in a manner similar to Katz, Weinberg & Hernquist (1996) (hence- 
forth KWH96). The code implements radiative cooling for an opti- 
cally thin gas of helium and hydrogen in collisional ionization equi- 
librium with a user-specified (time-dependent) ultraviolent back- 
ground. Radiative cooling allows the gas to dissipate energy and 
form dense, cold clumps, which as we describe below are prone 
to star formation. For the densities resolved in these simulations, 
radiative cooling due to two-body processes is important, namely, 
collisional excitation of neutral hydrogen (H°) and singly ionized 
helium (He + ), collisional ionization of H°, He and He + , standard 
recombination of H + , He + , and He ++ , dielectric recombination of 
He + and free-free emission. The adopted cooling rates are given 
in Table 1 of KWH96. The heating rate from photoionization can 
be computed once the intensity of the UV background is specified. 
Since this heating rate is very low inside disks in the local universe, 
we do not make use of this heating rate from a UV background. 

The version of GADGET-2 that we employ uses a sub- 
resolution model for energy injection from supernovae that is pro- 
portional to the star formation rate (Springel & Hernquist 2003). In 
this approach, the ISM is taken to be composed of cold clouds from 
which stars can form, which are embedded in a pressure-confining 
hot phase. The effective equation of state for this two-component 
gas (Eq. 24 in Springel et al. 2001) is taken to be 0.25 for our fidu- 
cial cases. In essence, employing such a equation of state for the 
star forming gas allows one to transition from the isothermal case 
to a stiffer than isothermal equation of state, where in the latter case 
energy injection from supernovae (which is proportional to the star 
formation rate) allows one to evolve gas-rich disks for many dy- 
namical times. We vary the equation of state in Appendix B and 
study its effect on the gas-star angular momentum transport. 

We discuss in this paper both simulations where star forma- 
tion is turned on and those where star formation is turned off. For 
simulations where star formation is turned on, we employ a density 
dependent Kennicutt-Schmitdt algorithm, to convert gas to stars in 
these simulations. The star formation rate is given by: 

where p is the total gas density, f3 ~ 0.1 for a Salpeter IMF, and 
the quantity t* is given by: 

/ \ -1/2 

**CP) = *o I — J . (3) 
\pthj 

where p t h is fixed by the multiphase model and is the threshold 
density above which the gas is assumed to be thermally unstable to 
the onset of a two-phase medium composed of cold star forming 
clouds embedded in a hot, pressure-confining phase (SH03). We 
take to, the gas consumption timescale, to be 4.5 Gyr and 8 Gyr 
which are intermediate values between that adopted by SH03 and 
Springel et al. (2005) and agree reasonably well (to within a factor 
of 2-3) with gas consumption time scales cited by Kennicutt (1998). 
We aim to highlight in this paper general trends in the dynamics 



that are seen across simulations. Due to the large dynamic range 
in global galaxy simulations, star formation is currently treated in 
a prescriptive way. To minimize dependence of our results on our 
choice of a star formation prescription, we also include results from 
simulations where the star formation is turned off. For the simula- 
tion where star formation is turned off, there is no energy injection 
from supernovae (as SH03's model takes the energy injection rate 
to be proportional to the star formation rate), and the gas remains 
isothermal at 10 4 K. We discuss in a future paper results from var- 
ious star formation prescriptions (where star formation is mediated 
by other factors such as shocks, magnetic fields and turbulence) for 
turning gas into stars. 

The initial conditions that we choose for the disks represent 
Milky Way type galaxies and are similar to those described in Cox 
et al. (2006) and Springel et al. (2005). Namely, we generate a 
stable isolated disk galaxy that is initially in equilibrium. The ex- 
tended dark matter halo is taken to have a Hernquist (1990) profile, 
with a concentration parameter of 10. We model the disk compo- 
nents of gas and stars as exponential disks, with the fiducial case 
having a scale length of 3 kpc. The primary simulation that we fo- 
cus on here has a virial velocity of 160 km/s, with a rotation curve 
and mass similar to the Milky Way (M tota i = 1.3 x 1O 12 M and a 
stellar mass of 5.4 x 1O 1O M0 initially). We do not include a bulge 
in the simulations reported here (there is no qualitative difference 
as far as our main results are concerned), and take the baryonic gas 
mass fraction to be 10 % for the fiducial case, with the gas dis- 
tributed initially in an axisymmetric exponential disk having the 
same scale length as the stars. (In retrospect, inclusion of an ex- 
tended gas disk as seen in observations by Wong & Blitz (2002) 
would have been more empirically suited) Values of these parame- 
ters for the other simulations are reported in Table 1 . These values 
yield an initial Toomre Q parameter for the stars that is slightly 
in excess of unity in the inner regions, and of order unity for the 
gas. We have performed resolution studies to confirm that simula- 
tions with the given number of particles and softening lengths yield 
converged results. The primary simulations that we analyze here 
are the collisionless SpNGNB simulation, its gaseous analogue 
(SpGNSFR), and its gaseous, star-forming analogues (SpGKS and 
SpGKS8). For a given cooling curve and star formation prescrip- 
tion, there are two aspects of the gas physics that can affect the 
gas-star dynamical coupling in these simulations, namely the ar- 
tificial bulk viscosity and the equation of state for the gas. These 
aspects of the simulation methodology and their effect on angular 
momentum transport and treated in Appendix A and B respectively. 



3 SIMULATIONS OF THE COLLISIONLESS 
COMPONENTS OF GALAXIES 

We depict in Figure 1 the projected (i.e., we take our stellar density 
in x,y,z and project it on the z-axis) stellar density (E) out to 2 Gyr 
for the SpNGNB simulation. This stellar disk is only somewhat 
susceptible to spiral structure as it is initialized with Q > 1 in the 
inner parts of the disk, with Q ~ 1 in the outer regions. A more 
quantitative representation of the spiral structure is given in Figure 
2, where we have plotted the m = 2 Fourier component of the 
stellar surface density. We obtain this by taking the amplitude of 
the real and imaginary components of a m =2{r), where a m (r, t) is 
given by: 

a m (r,t) = ± J*\(r,<b,t)e-^d0 . (4) 
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Table 1. Simulation Parameters 
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Figure 1. Projected Stellar Density images for the SpNGNB simulation out 
to 2 Gyr. The length of the box is 10 kpc. 
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Figure 3. The time evolution of J z for the stars (a) and the halo (b), for the 
SpNGNB simulation for < r < 7 kpc and 7 kpc < r < 13 kpc 



Figure 2. The m = 2 Fourier component of the projected stellar density at 
three time snaps for the SpNGNB simulation. 



Figure 2 shows that that the power in the m = 2 disturbance (rela- 
tive to the m=0 axisymmetric mode) varies by a factor of ~ 10 in 
magnitude peak-to-peak in the inner regions of the disk (r < 8 kpc) 
over the timescale of a Gyr. 

The spiral patterns disappear within an orbital time due to disk 
heating, an effect also seen by SC84. Rybicki (1971) estimated 
that the relaxation timescale is comparable to the orbital time for 
purely stellar razor-thin disks. The swing amplification mechanism 
(Toomre 1981) is effective for low values of Q - as the spiral pat- 
terns swing from leading to trailing in a shearing disk, they are 
amplified across co-rotation. In the absence of dissipation, this in- 
creases the velocity dispersion of the stars secularly until the disk 
is no longer susceptible to this amplification mechanism. 




.0 1.5 2.0 
time (Gyr) 

Figure 4. The time variation of the pattern speed, out to 3 Gyr (identified 
using Tremaine & Weinberg 1984) for the SpNGNB simulation (no gas, no 
bulge, stars and live DM halo only) 
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As noted earlier, the total angular momentum of the simula- 
tions is conserved. However, different galaxy components are free 
to exchange angular momenta. We show in Figure 3a the angular 
momentum (in arbitrary units) of the stars and dark matter halo 
particles within > r > 7 kpc and 7 kpc > r > 13 kpc (Fig- 
ure 3b). The stars inside co-rotation (which is at ~ 7 kpc) lose 
angular momentum to the halo; of particular note is the steady de- 
cline of angular momentum in the stellar disk. The massive slowly 
rotating halo is then a net sink of angular momentum, as has been 
found in previous simulations of bar instabilities with live halos 
(Athanassoula 2003; Martinez- Valpuesta et al. 2006; Berentzen et 
al. 2007). The basic idea that energy and angular momentum can 
be absorbed by the outer part of a galaxy if it is rotating slower 
than a given mode was discussed in the early studies of angular 
momentum transport in galaxies (Lynden-Bell & Kalnajs 1972). 

We show in Figure 4 the time variation of the pattern speed, 
computed following Tremaine & Weinberg (1984). This method 
derives the pattern speed on the assumption of a well-defined pat- 
tern speed and the tracer (E) obeying the continuity equation. This 
gives a pattern speed in terms of the projected surface densities 
(S(x, z)), which we have computed in the center of mass frame, 
and velocities (v(x, z)) of the stars: 

Q _ JZo H Z ) dz JZv E ( Z ' Z ) V ( X ' Z ) dX (5) 

S^ac, h(z)dz E(a;, z)xdx 

which we have expressed the pattern speed here for an edge- 
on disk; the weight function h(z)dz can be chosen arbitrarily. 
For this simulation, we find that the pattern speed is on average 
~ 27 km/s/kpc (for t>500 Myr), which places co-rotation at 
~ 7 kpc, and is in agreement with expectations from modal theory 
as we discuss later. The early-time (t = 250 Myr) pattern speed 
is high (which would imply a co-rotation radius close to the very 
inner regions of the disk), as the m = 2 spiral mode is not yet 
the dominant mode in most of the disk. For this collisionless sim- 
ulation, there is a 97 % change in the pattern speed out to 3 Gyr, 
and a 38 % change for t > 500 Myr. As we discuss later, simula- 
tions that include a gaseous component show less variance in their 
pattern speeds as a function of time. 



4 SELF-EXCITED QUASI-STEADY SPIRAL 

STRUCTURE - SIMULATIONS WITH DARK MATTER, 
STARS, AND GAS 

We now proceed to describing our results from simulations that 
include a gaseous component. We focus our analysis on angular 
momentum transport from the gas to the stars. The torque exerted 
by the self-gravitating gas disk on the stellar disk is given by - 

r.= r rs«(r,0 ^ r '*w , (6) 

Jo Jo d<t> 

where it is specifically the asymmetric components that contribute 
to this integral. 

We depict in Figure 5 the projected stellar density from the 
SpGNSFR simulation, which contains 10% of the disk mass in the 
gaseous component. Star formation is turned off in this simulation. 
As is clear, the stellar disk is unstable to spiral structure for a longer 
time when there is gas, and thin spidery asymmetric spiral arms 
are seen in the gaseous component (Figure 6). Figure 7 shows the 
time variation of the power in the m — 2 Fourier mode, and can be 
contrasted with Figure 2. In this case, when gas is included in the 
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Figure 5. Projected Stellar Density images for the SpGNSFR simulation 
out to 2 Gyr. The length of the box is 10 kpc. 




Figure 6. Projected Gas Density images for the SpGNSFR simulation out 
to 2 Gyr. The length of the box is 10 kpc. 



simulation, the m — 2 disturbance is roughly time-steady - there is 
less than a factor of three varation in the power peak-to-peak. 

Figure 8 depicts the evolution of angular momentum for the 
gas, stars, and halo interior and exterior to co-rotation. Angular mo- 
mentum transport from the gas to the stars interior to co-rotation is 
able to torque the stellar disk and serves to self-regulate the disk 
structure. Figure 8b shows the angular momentum of the stars is 
considerably more time-steady than its purely collisionless ana- 
logue (Figure 3a), as the stars are gaining the angular momentum 
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Figure 7. The m = 2 Fourier component of the projected stellar density at 
three time snaps for the SpGNSFR simulation. 
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Figure 9. Azimuthal phase shift in density (arbitrarily scaled) between the 
gas and the stars in the SpGNSFR simulation at 250 Myr at co-rotation 
(r = 7 kpc). 
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Figure 10. The time variation of the pattern speed, out to 2 Gyr (identified 
using Tremaine & Weinberg 1984) for the SpGNSFR simulation, with star 
formation turned off. 



o f 1 

12 3 

time (Gyr) 

Figure 8. The time evolution of J z for (a) the gas, (b) the stars; overplotted 
is J z for the stars for the SpNGNB simulation as a reference point for a case 
without gas (in this case the stars lose 7 % of their angular momentum inside 
of 7 kpc out to 2 Gyr in contrast to the 35 % loss in the purely collisionless 
case), and (c) the halo for the SpGNSFR simulation within < r < 7 kpc 
and 7 kpc < r < 13 kpc 



that is lost by the gas (Figure 8a). There continues to be a slight net 
transfer to the live halo. 

The sharp azimuthal response of the gas and the smooth re- 
sponse of the stars is shown in Figure 9 for the 250 Myr snapshot 
- this also shows that the gas leads the stars, with a phase shift of 
~ 0.5 radian close to co-rotation. Given that the potential traces 
the density, Figure 9 also shows why the gas is able to exert torques 
on the stars - ^$ gas S> ^$ s tars, as the gas is compressed and 
shocks, it forms the sharp azimuthal features shown in Figure 9. 
The response of the gas to a time-steady forcing potential in the 
nonlinear regime was studied by Shu, Milione & Roberts (1973) 
(henceforth SMR). SMR found that the response of the gas can be 
quite nonlinear even when the amplitude of the stellar spiral forc- 
ing is a small fraction of the axisymmetric field. The response of the 



gas is much sharper than that of the stars as the response of some 
galactic component to a forcing potential is proportional to the in- 
verse square of its dispersive speed. SMR also discussed in detail 
how the breadth of the zone of high gas compression depends on 
whether the Doppler shifted phase velocity of the stellar density 
wave (or the component of the circular motion that is perpendicu- 
lar to the wave front) is greater or less than the sound speed of the 
gaseous component, which results in narrow and broad zones of gas 
compression respectively. In this simulation, where star formation 
is turned off, there is no energy injection from supernovae. This 
allows the sound speed of the gas to be low enough to satisfy the 
criterion of SMR, and can be contrasted with the cases discussed 
later in this paper where star formation is turned on along with an 
injection of energy prescribed by the model of SH03 which results 
in broad gaseous arms, particularly in the outer parts of the galaxy. 
Time-dependent numerical calculations by CLS of the response of 
the gaseous disk when forced by a spiral stellar potential displayed 
similar sharp azimuthal features. 

Figure 10 displays the time variation of the pattern speed for 
the SpGNSFR simulation. The average pattern speed of this simu- 
lation is ~ 28 km/s/kpc, which places co-rotation at r ~ 7 kpc, 
and is in agreement with analytic expectations of the pattern speed 
from modal theory, i.e., a few scale lengths (Bertin et al. 1998) and 
from observational comparisons (Lin et al. 1969; Shu et al. 1970). 
There is no inner Lindblad resonance, and the outer Lindblad res- 
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Figure 11. Projected Stellar Density images for the SpGKS simulation out 
to 2 Gyr. The length of the box is 10 kpc. 




Figure 12. Projected Gas Density images for the SpGKS simulation out to 
2 Gyr. The length of the box is 10 kpc. 

onance is at ~ 15 kpc. We note that in contrast to the SpNGNB 
collisionless simulation, this simulation exhibits less variance in its 
pattern speed, with a maximum difference of 21 % out to ~ 2 Gyr. 

It is interesting to note that the average radial velocities of 
the gas (at times when the spiral structure is well established) are 
of order ~ 30% of the sound speed, which would cause gas to 
flow radially inwards ~ 2 kpc on the time scale of a Gyr. This 
effect cannot be neglected on the scale of the spiral arms, and the 
radial velocities in these simulations are in rough agreement with 
observational estimates (Wong et al. 2004). Radial inflow may also 
be produced in the outskirts of isolated galaxies due to the torque 
on the gas from new stars that exist along spiral arms far outside the 
optical radius as seen in recent Galex images (Thilker et al. 2007). 
We investigate the large scale transport of gas from extended HI 
disks as a possible mechanism for fueling star formation over many 
Gyr in a forthcoming paper. 

4.1 Star Formation Cases 

Figures 1 1 and 12 show the projected stellar and gas surface densi- 
ties for the SpGKS simulation, where star formation is prescribed 
via the Kennicutt-Schmidt relation with a gas consumption time 
scale of 4.5 Gyr. Although the results are qualitatively similar to 
the SpGNSFR case, the spiral structure in the gas is not as sharp 
due to energy injection from star formation as prescribed in the 
multiphase model of SH03 (which may be expected from the work 
of SMR as we have discussed in §4), and the stellar disk is also not 




Figure 13. Projected density for the new stars formed in the SpGKS sim- 
ulation out to 2 Gyr. The length of the box is 10 kpc. The first snapshot is 
blank as there is no star formation at t = 0. 
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Figure 14. The m = 2 Fourier component of the projected stellar density 
at three time snaps for the SpGKS simulation. 

as unstable to spiral structure over 2 Gyr. Figure 14 shows the time 
evolution of the m = 2 Fourier component of the stellar surface 
density, which is also in qualitative agreement with the SpGNSFR 
simulation, although the variation in power over a Gyr is somewhat 
larger. The time evolution of the angular momentum (Figures 15) of 
the stars is not as time steady as in the SpGNSFR simulation - over 
the period of 2 Gyr, the angular momentum of the stars declines by 
about 20 %, in contrast to the 7 % decrease for the SpGNSFR simu- 
lation, but both of these are below the 30% stellar angular momen- 
tum loss seen in the SpNGNB simulation that did not include gas. 
Of note is the projected density of the new stars as shown in Figure 
13, which is more akin to the sharper response of the gas shown 
in Figure 12. Figure 16 depicts the density vs azimuth at 7 kpc for 
the SpGKS simulation and is similar to the SpGNSFR case, though 
with a slightly broader azimuthal response in the gaseous compo- 
nent. The time variation of the pattern speed is shown in Figure 17, 
with an average pattern speed of ~ 28 km/s/kpc, which places 
co-rotation at ~ 7 kpc. There is a maximum variance of 27 % in 
the pattern speed for the SpGKS simulation out to 3 Gyr. 

Gas consumption time scales in nearby spiral galaxies are in- 
ferred to be on the order of a few Gyr (Kennicutt 1983; Wong & 
Blitz 2002; Crosthwaite & Turner 2007), with a range between 1- 
10 Gyr (Wong & Blitz 2002). The values are uncertain both due to 
use of various tracers to infer the rate of star formation and hence 
the gas depletion time scales, as well as due to the necessary use 
of globally averaged measures which recent works have been able 
to address observationally (see e.g., Kennicutt et al. 2007). Our use 
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Figure 15. The time evolution of J z for (a) the gas, (b) the stars shown 
along with the comparison to the purely collisionless SpNGNB case (in this 
case the stars lose 38 % of their angular momentum inside of 7 kpc out 
to 3 Gyr in contrast to the 43 % loss seen in the purely collisionless case), 
and (c) the halo for the SpGKS (in my notation SpGKS) simulation within 
< r < 7 kpc and 7 kpc < r < 13 kpc 
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Figure 16. Azimuthal phase shift in density (arbitrarily scaled) between the 
gas and the stars in the SpGKS simulation at co-rotation at 250 Myr. 



100 



- 1 — 1 — • — • — • — I - 



1 I ^_ 

0.0 



_l . , , , 1_ 



0.5 



1.0 1.5 2.0 
time (Gyr) 



2.5 



3.0 



Figure 17. The time variation of the pattern speed, out to 3 Gyr (identified 
using Tremaine & Weinberg 1 984) for the SpGKS simulation. 



what longer (by ~ 2) gas consumption time scale for this simu- 
lation allows the gaseous component to couple dynamically to the 
stellar disk for a slightly longer time. There is likely to be con- 
tinued accretion of gas from the intergalactic medium, which may 
serve to replenish the spiral instability over even longer timescales, 
a point that we discuss qualitatively in §5. Finally, Figure 24 display 
the time variation of the pattern speed for the SpGKS8 simulation, 
which is on average 27 km/s/kpc and varies by ~ 40% over 3 



Gyr. 



of a density dependent star formation prescription is simplistic; star 
formation in real galaxies may be mediated magnetically (Shu et al. 
2007) and or through turbulence (Krumholz & McKee 2005), and 
may be understood in a simple way in terms of the pressure (Blitz 
& Rosolowsky 2006). Nonetheless, use of this standard star for- 
mation prescription is fairly pervasive (SH03; Springel et al. 2005; 
Cox et al. 2006) and does agree with globally averaged measures 
of the star formation rate in a diverse range of galaxy types. 

We also include the SpGKS8 case, with a gas consumption 
time scale of 8 Gyr. Figures 18-20 display the projected images 
stellar, gas, and new stars respectively, which are similar to the 
SpGKS simulation, but with a slightly greater spiral disturbance 
seen in Figure 18 at t = 1.5 Gyr, which is also reflected in the late 
time m — 2 Fourier component shown in Figure 21. The some- 
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Figure 18. Projected Stellar Density images for the SpGKS 8 simulation out 
to 2 Gyr. The length of the box is 10 kpc. 
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Figure 19. Projected Gas Density images for the SpGKS8 simulation out to 
2 Gyr. The length of the box is 10 kpc. 




Figure 20. Projected density for the new stars formed in the SpGKS8 sim- 
ulation out to 2 Gyr. The length of the box is 10 kpc. The first snapshot is 
blank as there is no star formation at t = 0. 



5 DISCUSSION 

A number of papers have discussed the role of gas in providing a 
dissipative and therefore a saturation mechanism (Roberts & Shu 
1972; Lubow, Balbus & Cowie 1986; Bertin et al. 1998). Our time- 
dependent simulations treat the gas-star interplay self-consistently. 
The gaseous component is able to cool radiatively and dissipate 
energy, and thereby remain cool enough to be prone to the spiral 
instability; the azimuthal phase shift between the gas and the stars 
allows for the donation of angular momentum from the gas to the 
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Figure 22. The time evolution of J z for (a) the gas, (b) the stars overplotted 
with the reference SpNGNB case (in this case the stars lose 23 % of their 
angular momentum inside of 7 kpc out to 3 Gyr in contrast to the 43 % 
loss seen in the purely collisionless case), and (c) the halo for the SpGKS8 
simulation within < r < 7 kpc and 7 kpc < r < 13 kpc 
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Figure 21. The m = 2 Fourier component of the projected stellar density 
at three time snaps for the SpGKS8 simulation. 



stellar disk. We find that this expropriation of angular momentum is 
a critical factor in the long-term persistence of spiral structure in the 
stellar disks of galaxies. Specifically, the torques that the gaseous 
spiral sub-structure (those features that are not azimuthally sym- 
metric) exert on the stellar disk transfer angular momentum and 
therefore lead to a roughly time-steady spiral structure in the stel- 
lar disk. The inclusion of a cold gaseous component also reduces 
the effective Q parameter of the two-fluid system (Jog & Solomon 
1984; Bertin et al. 1998; CLS03). Springel et al. (2005) elucidated 
the effects of the Springel & Hernquist (2003) multiphase model 
by performing a parameter space survey; due to the energy injec- 
tion from supernovae that is implicit in this model, extremely high 
gas fractions are not unstable to spiral structure. We discuss this 
point further in Appendix B. 

Kendall et al. (2008) recently produced a composite surface 
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Figure 23. Azimuthal phase shift in density (arbitrarily scaled) between the 
gas and the stars in the SpGKS8 simulation at co-rotation at 250 Myr. 
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Figure 24. The time variation of the pattern speed, out to 3 Gyr (identified 
using Tremaine & Weinberg 1984) for the SpGKS8 simulation 



density map of M81 using Spitzer's IRAC 3.6 /im, 4.5 fim, 8 fim, 
optical, B, V, I and 2MASS K B band data. They extract the am- 
plitude and phase of the density from the map, and find that the 
response of the gas, which is traced by dust emission in the 8 /im 
IRAC band, is phase-shifted in azimuth relative to the stellar den- 
sity wave. Our analysis of the simulations here is consistent with 
their findings, with azimuthal offsets between the gas and stars of 
order one radian. In a future paper, we perform radiative transfer 
calculations through the simulations to calculate simulated images 
as done previously by Chakrabarti & Whitney (2008) for ULIRGs 
in the continuum bands. Our gas surface density images (Figure 
6 for example) show a much sharper response as seen in B-band 
images of spirals, while the stellar surface density images (Figure 5 
for example) show a much smoother response as seen in the K band 
images of local spirals (Block et al. 1996; Elmegreen & Elmegreen 
1984; Rix & Zaritszky 1995). 

A major caveat to our simulations here is that we have not in- 
cluded the cosmological infall of gas, which can potentially serve to 
replenish the supply of gas over Gyr timescales as the pre-existing 
supply of gas is steadily exhausted by star formation and trans- 
ported inwards by gravitational torques. Colavitti et al. (2008) have 
recently analyzed the chemical evolution of a spiral galaxy in a cos- 
mological context and find a gas infall prescription that reproduces 
observed chemical properties. Such a prescription may also be used 
to study the long term dynamics of spiral structure in a more cos- 
mologically motivated way than we have done here. As we have 
studied galaxies evolving in isolation which cannot make use of the 



continued accretion of gas to replenish the critical supply of a cold 
dissipative component, as gas is tranported inwards or consumed 
by star formation, we very likely find a lower limit (~ 2 Gyr) for 
the lifetime of spiral structure in disk galaxies. This time scale may 
well be extended by fresh infall of gas. Gas-rich mergers will cer- 
tainly affect the stability of disks; however, detailed knowledge of 
their frequency, mass ratio (i.e., minor or major mergers), and fi- 
nal dynamical impact cannot be obtained without performing high 
resolution hydrodynamical cosmological simulations. Recent ob- 
servational works estimate that a third of LIRGs (Lir > lO 11 !/©) 
are produced by major mergers at z ~ 1, while the rest are mor- 
phologically akin to spirals (Elbaz et al. 2007). The morphological 
variation of normal galaxies with redshift (due to the influence of 
other factors such as mergers) is another issue that we have not 
treated here. Our approach here has been to study the simple case 
of a galaxy evolving in isolation and to understand specifically the 
gas-star dynamical coupling that is a critical factor in the long term 
persistence of spiral structure in stellar disks. Due to the relative 
simplicity of our approach, we have been able to study aspects of 
the gas physics that impact the angular momentum transfer from 
the gas to the stars, such as the effect of star formation (by studying 
simulations where star formation is turned off, i.e., in §4, as well 
as simulations where star formation is turned on and the gas con- 
sumption timescale is varied, i.e., in §4.1), as well as the effect of 
the artificial bulk viscosity (Appendix A), and the equation of state 
(Appendix B). 



6 CONCLUSION 

• We find that angular momentum exchange between the self- 
gravitating gas disk and the stars leads to relatively long-lived spiral 
structure in the stellar disk, with spiral structure lifetimes of order a 
few Gyr. Spiral structure is washed out, when it is incipient, in sim- 
ulations that do not include a cold, dissipative gaseous component. 

• The gas disk is able torque the stellar disk and donate an- 
gular momentum interior to co-rotation. The gas is able to cool ra- 
diatively and dissipate energy - leading to a time lag (or azimuthal 
phase shift on the order of ~ 1 radian close to co-rotation) between 
the response of the gas and the response of the stars. Even though 
the gas mass is less than the stellar mass, the gas can exert appre- 
ciable torques (and donate non-negligible angular momentum) as 
the azimuthal response which is highly asymmetric (or azimuthal 
gradient of the potential) is much larger than that of the stars. Over 
timescales of many Gyr, there is radial inflow of gas which ulti- 
mately limits the donation from the gas (if it is not replenished 
due to cosmological infall); star formation (depending on gas con- 
sumption time-scales) will convert the cold, dissipative component 
to stars which heat up and have no cooling vent. 

• The primary (specifically the m = 2) Fourier components 
of the stellar surface density are roughly time-steady in simulations 
that include a gaseous component, with variation in surface density 
of ~ 3 over a Gyr. The pattern speed of the m — 2 component is 
more time steady in simulations that include gas than those that do 
not. 

• The halo is a net sink of angular momentum, as has been 
seen in previous studies of bar instabilities with live halos. 
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APPENDIX A: ANGULAR MOMENTUM TRANSFER & 
ARTIFICIAL BULK VISCOSITY 

We discuss here the consequence of varying the artificial bulk vis- 
cosity parameter and its effect on angular momentum transfer from 
the gas to the stars. Figure 23 (a-c) depicts the angular momen- 
tum transfer for the SpGNSFR simulation with the bulk viscosity 
parameter set to zero. This is similar to the collisionless case, and 
results in the angular momentum of the stars steadily declining, 
as there is no transfer from the gas due to a lack of phase shift 
between the gaseous and stellar components. The standard value 
for the artificial bulk viscosity that we have used in the paper is 
0.75; this case is shown in Figures 24 (a-c). Figures 25 (a-c) depict 
the SpGNSFR case with a artificial bulk viscosity of unity, which 
results in a very similar behavior of angular momentum transfer 
as our standard case. The zero viscosity simulation (SpGNSFRO) 
shows a 24 % loss in the stellar angular momentum out to 1.4 Gyr 
inside of 7 kpc (which is comparable to the 30 % loss in stellar 
angular momentum for the reference SpNGNB purely collision- 
less case), while the standard viscosity case (SpGNSFR) and the 
SpGNSFRl both experience a 6 % loss out to 1.4 Gyr interior to 7 
kpc. In summary, the effect of lowering the artificial bulk viscosity 
is to approximate the purely collisionless problem, thereby remov- 
ing the phase shift in azimuth that is needed for the gas to exert a 
torque on the stellar component. Aside from the effect of viscosity, 
gas particles also exert pressure - this will not affect the angular 
momentum transport directly, but will affect the propensity of a gas 
disk to be prone to spiral structure. We discuss the effect of pressure 
(i.e., assumptions about the equation of state) in Appendix B. 
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Figure Al. J transfer for SpGNSFR simulation with the bulk artificial vis- 
cosity set to zero (a) gas (b) stars, overplotted with the reference SpNGNB 
case which has a 30 % loss in angular momentum out to 1 .4 Gyr interior to 
7 kpc compared to the 24 % loss for this zero viscosity case, and (c) halo 
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Figure A2. The time evolution of J z for (a) the gas, (b) the stars over- 
plotted with the reference SpNGNB case which has a 30 % loss in angular 
momentum in contrast to the 6 % loss seen in this standard viscosity case, 
and (c) the halo for the SpGNSFR simulation within < r < 7 kpc and 
7 kpc < r < 13 kpc. This is the standard case with the bulk artificial 
viscosity set to 0.75. 



Figure A3. J transfer for SpGNSFR simulation with the bulk artificial vis- 
cosity set to one (a) gas (b) stars overplotted with the reference SpNGNB 
case which has a 30 % loss out to 1.4 Gyr in contrast to the 6 % loss in 
angular momentum seen here, and (c) halo 
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APPENDIX B: EQUATION OF STATE 

Another aspect of the gas physics that affects the stability of the 
gas disk, and indirectly impacts the angular momentum transport 
to the stellar disk is the assumed equation of state for the gas. As 
discussed by Springel et al. (2005), the equation of state, when it 
is taken to be greater than zero which is the isothermal limit, mod- 
els the effect of supernova feedback by over-pressurizing regions 
relative to an isothermal equation of state. In the isothermal limit, 
the pressure in compressed regions increases in proportion to the 
density. The multiphase model mediates that for p > p t h, the gas 
is thermally unstable and stars can form; the energy injection rate 
from supernovae in this model is proportional to the star formation 
rate. The compressed gas then would have an effective equation 
of state that is stiffer than isothermal if it retains this feedback en- 
ergy. One can expect that simulations performed with equations of 
state that are stiffer than isothermal are somewhat less prone to a 
spiral instability as the effective sound speed has been increased. 
We vary the equation of state from 0-1 in Figures 26(a-c) and find 
that the simulation performed with the isothermal equation of state 
(EQSOLR) shows the lowest decrease in stellar angular momen- 
tum, while the full multiphase model (EQS1LR) has a 25 % de- 
crease in the stellar angular momentum out to 1.4 Gyr. However, 
this relative decease - 18 % versus 25 %, as the equation of state 
varies from isothermal to the fully pressurized multiphase model 
is a small enough change to suggest that assumptions about the 
equation of state for the gas only slightly affect the angular mo- 
mentum transport to the stellar disk, for a gas fraction of 10 %. 
Springel et al. (2005) have performed a comprehensive parameter 
study in terms of the gas fraction and the equation of state to study 
the effects these parameters have on the gaseous disk. They find 
that for the full multiphase model (stiff equation of state) even pure 
gas disks are stable to fragmentation, with an increase in gas frac- 
tion lowering the susceptibility of the gas disk to spiral structure. 
In the multiphase model, increasing the gas fraction increases the 
star formation (as mediated by the adopted Kennicutt-Schmidt al- 
gorithm) and thereby the supernova feedback energy (which is pro- 
portional to the star formation rate) that is deposited into the gas. 
This can cause the effective temperature of the gas to be dynam- 
ically hotter than the stars. On the other exterme, extremely high 
gas fraction disks that are evolved with an isothermal equation of 
state undergo violent fragmentation. Thus, extremely high gas frac- 
tion disks (whether they are evolved with an isothermal equation of 
state or full multiphase model) will not be effective as angular mo- 
mentum donors to the stellar disk. 
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Figure Bl. Stellar angular momentum as a function of time for (a) the 
EQSOLR simulation with the equation of state to zero (isothermal), which 
shows a 18 % loss out to 1.4 Gyr (b) for the EQSSTLR simulation with the 
equation of state set to 0.25 (our standard case which is over-pressurized 
relative to the isothermal case); this has a 22 % loss in the stellar angular 
momentum out to 1.4 Gyr (c) for the EQS1LR simulation with the equation 
of state set to one that is the full multi-phase model shows a 25 % loss in 
angular momentum out to 1 .4 Gyr 
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